Mode-coupling theory and molecular dynamics simulation for heat conduction in a 
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We study heat conduction in a one- dimensional chain of particles with longitudinal as well as 
transverse motions. The particles are connected by two-dimensional harmonic springs together 
with bending angle interactions. The problem is analyzed by mode-coupling theory and compared 
with molecular dynamics. We find very good, quantitative agreement for the damping of modes 
between a full mode-coupling theory and molecular dynamics result, and a simplified mode-coupling 
theory gives qualitative description of the damping. The theories predict generically that thermal 
conductance diverges as iV^/^ as the size A'^ increases for systems terminated with heat baths at the 
ends. The A'"^''^ dependence is also observed in molecular dynamics which we attribute to crossover 
effect. 

PACS numbers: 44.10.-|-i, 05.45.-a, 05.70.Ln, 66.70.-t-f 



I. INTRODUCTION 

The problem of heat conduction is a well-studied field. 
More than two centuries ago, Joseph Fourier summarized 
the behavior of heat conduction by the law that bears 
his name. This law describes phenomenologically that 
the heat current is proportional to the temperature gra- 
dient. The detailed atomistic theories of heat conduc- 
tion appeared only much later. For heat conduction in 
gas, the simple kinetic theory gives the result k = ^Cvl, 
where C is specific heat, v is average velocity and I is 
mean-free path. Peierls' theory of heat conduction in in- 
sulating solids j3] is a classic on this subject. These early 
theories deal with mostly the relevant three dimensions. 
It turns out that low dimensional systems are more in- 
teresting and in some sense strange. An analysis of a 
simple one-dimensional (ID) harmonic oscillator model 
shows that there is no well-defined temperature gradi- 
ent, the thermal conductivity diverges with system sizes 
as N^^'^ or N, depending on the boundary conditions. 
There is a general argument, for momentum conserving 
systems, that the thermal conduction in ID is necessarily 
divergent 

There have been many analytical and numerical stud- 
ies of ID heat conduction (see ref. ^ and for review). 
We'll mention some of the most relevant papers to cur- 
rent work. The work of Lepri et al by mode-coupling 
theory and molecular dynamics suggests a divergent ther- 
mal conductivity exponent of 2/5, i.e., k cx 

for a ID 

chain model with Fermi-Pasta-Ulam (FPU) interactions. 
Mode-coupling theory is usually applied in the dynam- 
ics of liquids 1^, Q . The first use of this theory in the 
context of heat conduction appears only recently, mostly 
due to Lepri and his collaborators JJ . Perverzev analyzed 
the same problem with the Peierls theory for phonon gas 
and gave the same conclusion of 2/5 exponent The 



result of 2/5 is also supported bynumerical simulation 
from several groups [H m m These results are 
supposed to be universal to some extent. However, it is 
challenged by a different result of 1 /3 by Narayan and Ra- 
maswamy based on fluctuating hydrodynamics and 
renormalization group analysis. The numerical result for 
this 1/3 law is not convincing, as for the same model - 
the hard-particle gas model - some obtained 1/3 ^16., 17|. 
while others obtained different value 1/4 |0| . But for the 
FPU model, there is no good evidence for an exponent 
of 1/3 19]. 

When momentum conservation is broken such as the 
one with on-site potential, the heat conduction can be- 
come normal again like the Frenkcl-Kontorova model 
and the 0* model [11 ij. 

In order to understand the underlying microscopic dy- 
namical mechanism of the Fourier law, a different class 
of models - billiard channels - have been introduced and 
studied in recent years [2^ |23 | . Various exponent values 
are found in such systems. Thus, it is believed that a 
universal constant does not exist at all. Instead, the di- 
vergent (convergent) exponent of the thermal conductiv- 
ity is found related to the power of super (sub) diffusion 



Besides the theoretical significance of heat conduction 
research in low dimensional systems, it is also of practical 
importance. Recent development of nanotechnology will 
enable us to manufacture devices with feature sizes at 
molecular level. The understanding of heat conduction 
mechanism will allow us to control and manipulate heat 
current, and eventually to design novel thermal devices 
with certain function p3 |. To this end, more realistic 
physical models are necessary. Among many others, nan- 
otubes and polymer chains are most promising. There 
have been a number of numerical works to com pute the 
thermal conductivity of the Carbon nanotubes |25l |2^ . 
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Recent molecular dynamics (MD) study of Carbon nan- 
otubes with realistic interaction potential suggested a di- 
vergent thermal conductivity for narrow diameter tubes 
[^T", '28] . The quantum effect of such systems is also very 
interesting _29, 30] . 

We study the heat conduction of a ID solid, as a clas- 
sical system. A brief version of this paper is reported 
in ref. [IJ. In the rest of the paper, we introduce the 
quasi-one-dimensional chain model in Sec. |nj We dis- 
cuss the basis of the mode-coupling theory, the projection 
method in Sec. IIIII The mode-coupling approximations 
and their numerical and analytical solutions are discussed 
in Sec. II VI and Ivl The basic output of the mode-couphng 
analysis is the dependence of damping of the modes with 
the wave-vector of the modes. We find that the trans- 
verse modes are diffusive, with (x p^, while the lon- 
gitudinal modes are super-diffusive, 7]! oc p'^/^, where 7^ 
is the decay rate for mode with momentum or lattice 
wave number p. We discuss the relationship between the 
damping of the modes with the heat conductivity through 
Green-Kubo formula in Sec. I VII Our mode-coupling the- 
ory predicts that the heat conductance diverges with the 
1/3 exponent when the transverse motions are impor- 
tant, while 2/5 is recovered if the transverse motions can 
be neglected. In Sec. IVIII we present nonequilibrium 
molecular dynamics results (with heat baths) of heat con- 
ductance and compare with mode-coupling theory. We 
conclude in the last section. 



II. CHAIN MODEL 

Most of the previous studies considered only strictly 
ID models, with the Fermi-Pasta-Ulam (FPU) model as 
the most representative. The strictly ID models may not 
be applicable to real systems such as the nanotubes. Real 
systems of nanotubes or wires live in three-dimensional 
space. The added transverse motion and the flexibility 
of the tube at long length scales will certainly scatter 
phonons, and thus should have a profound effect on ther- 
mal transport. 

While a direct simulation of a realistic system such as a 
polyethelene chain with empirical force fields such as that 
in refs. [3^ l33| or nanotube with Tersoff potential [s^l 
is possible, we think it is useful to consider a simplified 
model which captures one of the important features of the 
real systems - transverse degrees of freedom. Therefore 
we propose to study the following chain model in two 
dimensions l31l: 



i i 

7^0^ COS 4, 



(1) 



where the position vector r = {x,y) and momentum vec- 
tor p = {px,Py) are two-dimensional; a is lattice con- 
stant. The minimum energy state is at (ia,0) for i = 



to — 1. If the system is restricted to yi = (corre- 
sponding to — 4-00), it is essentially a ID gas with 
harmonic interaction. The coupling Kr is the spring con- 
stant; signifies bending or flexibility of the chain, 
while (f)i is the bond angle formed with two neighboring 
sites, cos0i = — rii-i-rii, and unit vector — Ari/|Ari|, 
Ar^ = r.,+1 - r.,. 

Unlike the FPU model which does not have an energy 
scale, the second bond-angle bending term introduces an 
energy scale. In this work, we take mass m = 1, spring 
constant Kr ~ 1, and the Boltzmann constant fc^ = 1, 
thus the most important parameters are AT^ and temper- 
ature T. 



III. PROJECTION METHOD 

A. Basic theory of projection 

We follow the formulation of the projection method in 
ref. m. Let 



A 



0.2 

\anj 



(2) 



be a column vector of n components of some arbitrary 
functions of dynamical variables {p, q) . Each of the func- 
tion aj{p,q) can be complex. Later, we shall choose aj 
to be the canonical coordinates of the system. We use 
A^ — (a|, flj, • • ■ , a* ) to denote the Hermitian conjugate 
of A. The equation of motion for A is 



At = CAt 



aj{t) = Caj{t), 



where C is the Liouvillc operator 

^dH d ^dH d 
£ = - > h > . 

— ^ Oq Op — ^ 



dp dq 



(3) 



(4) 



Equation Q can be viewed as a partial differential equa- 
tion with variables {p,q) and time t. At = A(j)t,qt) = 
A(t^p,q), i.e., the quantity A at time t when the initial 
condition at t = is (p, q). Quantity without a subscript 
t will be understood to be evaluated at time t = 0, e.g., 
p — pt=a — p(0). A formal solution to Eq. Q is simply 
At = e^^A{p, q). 

The projection operator on a column vector X is de- 
fined by 



VX ^ {X,A^){A,A^)'^A, 



(5) 



where {X, A^) and {A, A^) are nxn matrices. The angu- 
lar brackets denote thermodynamical average in a canon- 
ical ensemble at temperature T. The comma separating 
the two terms is immaterial, but we use a notation of in- 
ner product. One can verify that V is indeed a projection 
operator, i.e., = V. 



3 



If we apply the projection operator V and V' = 1 — V 
to the equation of motion, we get two coupled equations. 
Solving formally the second equation associated with 7-", 
and substituting it back into the first equation, we ob- 
tain an equation for the projected variable that formally 
resembles a Langevin equation: 



T{t - s)As ds + Rt 



where i = \/"-T is the complex unit, and 

r(t) - {RuRl){A,A^)-\ 
= (i,At)(A,At>-i, 
i?t = e*^'^i?o, Ro = A-iQA = V'CA. 



(6) 



(7) 
(8) 
(9) 



This set of equations is deterministic and formally exact. 
The only assumption made is that equilibrium distribu- 
tion can be realized. 

What is more important is the correlation functions, 
which are physical observables. We define the normalized 
correlation function (correlation matrix) as 



Git)^{A„Al){A,A^)-\ 



(10) 



It is an identity matrix at i = and has the property 
that G{0) = iQ. From Eq. © we have. 



G{t) = inG{t) 



Tit - s)G{s)ds. 



(11) 



This equation can be solved formally using proper initial 
condition with a Fourier-Laplace transform. 



G[z] = r 

Jo 



'G{t)dt, 



and similarly for T(t). The solution is 

G[z] = {i{z-n)+T[z]y\ 



To simplify notation, we have used parentheses (e.g., 
G{t)) for functions in time domain, and square brackets 
with the same symbol (e.g., G[z]) for their corresponding 
Fourier-Laplace transform in frequency domain. 

The information about the system is in the memory 
kernel T{t). However, such a correlation function is dif- 
ficult to calculate, since the evolution of the "random 
force" Rt does not follow the dynamics of the original 
Hamiltonian system. For example, it is impossible to 
compute directly Rt from molecular dynamics. For this 
reason, a "true force" can be introduced which obeys the 
normal evolution, i.e., 



Ft=e'^Ro=At-inAt. 
The correlation function of the true force, 
TFit) = {Ft^F^){A,A^)-\ 



(14) 



(15) 



and that of the random force are related in Fourier space 
as |35i 



T[z]- 



Tpizr^ - {t{z-n)y 



(16) 



This completes the formal theory of projection due orig- 
inally to Zwanzig [s^ and Mori !37j. These results are 
formal and exact. They give us relation between correla- 
tion of the "force" and correlation of dynamical variables. 
They are the starting point for mode-coupling theory. In 
the next subsections, we apply it to our chain model and 
introduce a series of approximations to solve it. 



B. The chain model 

We now apply the projection method to our chain 
model. We choose the normal modes as the basic quanti- 
ties Uj that we are going to project out. There are several 
reasons for choosing the normal modes. To zeroth order 
approximation, each mode is nearly independent. The 
slowest process corresponds to long wave-length modes. 
The effect of short wave-length modes can be treated as 
stochastic noise (the random force Rt). 

We choose A to be the complete set of canonical mo- 
menta and coordinates: 



, A: = 0,1,--- ,7V-1, 



(17) 



(12) where 



(13) 



Ql 
Qi 



N-l 

m X ^ 



JV-l 



J=0 

N-l 



pt = Qi = y\ 



Pj,v' 



,i2Trjk/N 



J2Trjk/N 



(18) 
(19) 
(20) 
(21) 



3=0 



We have defined the position vector rj — [xj ,yj) — 
[uj + aj, Uj), so that uj and j/j are deviations from zero- 
temperature equilibrium position. Because the Fourier 
transform is a periodic function, the index k is unique 
only modulo N. As a result, we can also let k vary in the 
range -N/2 <k< N/2. We also note that Ql = Q-k. 

With these definitions, we can compute the matrix 
and expression for the true force F in the general theory. 
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We find that for {A, A^), the components are 



{pi:Qk'*) = 0, 



/3 = 



(22) 
(23) 
(24) 



We have used equal-partition theorem for the average 
kinetic energy expression, and the last equation merely 
defines the effective frequencies a)^ for each mode. Note 
that due to translational invariance, correlation between 
different k modes vanishes. Correlation between trans- 
verse and longitudinal modes also vanishes due to the 
reflection symmetry of Uj ^yj for the chain. Thus 
equal-time correlation for A is diagonal, 



(A At) 



1 // 

Cj-^ 



(25) 



We have defined (D as a 2N x 2N diagonal matrix with 



elements w^; / is a 2N x 2N identity matrix. Similarly, 
the correlation {A,A^) is found from 

0, 



{Q'kPk'l 



Shk'Si, 



1 

1 

7' 



(26) 
(27) 

(28) 
(29) 



The second equation is from a general virial theorem [38 
We have 







(30) 



The expression for the true force F ^ A — iilA is then 

(31) 



F 



Pi: 







Note that only the momentum sector has a nonzero value, 
and = is given approximately by Eq. (|50l) and 
(|51|l below. With this special form of F, the damping 
matrix T is also diagonal and is nonzero only in the PP 
components. With these results, Eq. H13|) becomes 



G[z] = 



izd 
d [iz 



l-f[z])rf 



(32) 



where d = (cj^ - z'^I + izf[z\)~'^ is a 27V x 27V diag- 
onal matrix, and T[z\ is the Fourier-Laplace transform 
of the correlation P{R[ , {R^)^), which is also diagonal 
(we'll denote as r^(t)). In particular, we have the usual 
expression for the normalized coordinate correlation. 



UJ 



z^ + 



zzT^iz] 



(33) 
(34) 



For simplicity, we drop the QQ subscript for the coordi- 
nate correlation for the rest of the paper. 

Finally, the relation between the random force corre- 
lation and true force correlation, Eq. (|16|) . becomes, 



1 



1 



(35) 



The force correlation j, can be computed from the 
correlation function P{Ff , (F^Y) . It is convenient to 
separate the linear term in the force from the nonlinear 
contribution. So we write 

F^it) = u^Qit) + Q{t) = (^2 _ ^l)Qit) + j^it)^ (36) 

where Cj is effective angular frequency and wq is bare 
angular frequency of the mode. We have dropped the 
indices k and ^ since these equations apply for any of the 
modes, /at is at least quadratic in Q. We note that the 
correlation of F^ [t) is linearly related to the coordinate 
correlation function 



Tf[z] 



-IZ{&^ - z2) (c^2 _ ^2)2 



(37) 



This is simply a consequence of Eqs. (I34|l and (|35|l . but 
can also be derived directly from the definition. The 
second derivative of Q(t) in frequency domain is Q[z] 
multiplied by (iz)'^. Using the fact that 



lim / Q(t)e~ 
«--o+ Jo 



'^^dt = ~z^Q[z] - Q(0) - izQ{0), 

(38) 

we can also derive Eq. 1)3 7() with the understanding that 
(• • • ) is an average over the initial conditions. Since g[z] 
is finite or at least should not diverge precisely a,t z = uj, 
this implies Fp-p] = 0. 

With the above results, we can derive an expression 
of the true force correlation in terms of nonlinear force 
correlation. 



= (F^(t)F^(O)*) 

= ACj^{Q{t)Q*iO)) + {fr,it)fN{Or) 



where Alj^ = loq — uP'. The mixed term can be expressed 
in terms of g[z\ by noting that /at = Q -\- tUgQ. The two 
mixed terms {fN{t)Q* (0)) and {Qit)f^{0)) are equal due 
to time-reversal symmetry. We find 

/"OO 

pCu^ / {fN{t)Q*me''''dt = {col - z^)g[z] - iz. (40) 



Finally, we have 



Tp[z] 



((2^2 _ ^2 _ ^2)^[^j ^ 2iz) + Tn[z]. (41) 
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We can also express g[z] in terms of the nonlinear part 
of the force correlation, Tjv (= /3(/Ar(i)/Ar(0)*)): 



iz{2LLiQ 



2\2 



(42) 



Again, since g[z\ cannot diverge precisely at z = loq^ this 
implies that F^vf^] must take a special form to cancel 
the apparent divergence. Thus, if we do not take care of 
these superficial divergences, we will not be able to make 
correct prediction for the correlation function. 

We can also relate the original damping function F to 
the nonlinear one, 



FM = 



oJq — HP-z"^ — izilP-V]S!\z\ 



(43) 



when taking the derivatives. A slightly modified Hamil- 
ton's equation (because of the use of complex numbers) 
describes the dynamics: 



dH 
OH 



This gives the following equations of motion: 



Qk = -^fc <?! 



where 



C^'^" Qk' Qk" ; 



Qi 



k'+k" = k 



■'k^Qk + ^ ^k'^k"Qt'Qk"^ 



k'+k"=k 



(48) 
(49) 

(50) 
(51) 



This last equation is useful for approximating the damp- 
ing function. All of these relations are exact. This ends 
our formal application of the projection method to the 
chain model. 



IV. MODE-COUPLING THEORY 

To make some progress for analytic and numerical 
treatment, we have to make some approximations. First, 
we'll consider small oscillations valid at relatively low 
temperatures. An approximate Hamiltonian for small 
oscillations near zero-temperature equilibrium position, 
keeping only leading cubic non-linearity in the Hamilto- 
nian, is then given by 



HiP,Q) - \ E {Pi:P'-k + <'QtQ'-k) 

k\^—\\ ,_L 

+ Ck^p^qQ^Qp Qq , 

k-\-p+q=0 mod TV 



where 



Il2 



sm , 

m N 



.-L2 



(44) 

(45) 
(46) 



are the 'bare' dispersion relations, and 
8« . k-K 



^k,p,q 



TTK sm — sm — sm — —a K, 

1/2 N N N\2 



+K^{-2 - 



27rp 2T:q 
.cos— cos— ) 



(47) 



The absence of Q"(3"(5'I term in Hamiltonian H44I) is due 
to the quadratic nature of the potential, while the ab- 
sence of the terms of the form Q^Q-^Q-^ and Q^^Q^^Q-^ is 
due to the reflective symmetry about y-axis of the Hamil- 
tonian. We view k and —A: as independent component 



^k,p 



-■k,p 



= ii 
Kr 



1 , kn , pn , (k + p)tt 
sm — sm — sm 

Nrm'a^ N N N 

2 , 2kT: 2p7r. 

^X4-2-Hcos— -t-cos— ^ 



(52) 



1 



Af-2 + 



sm — sm - — sm ^ — — I Kr 

" N N \ 



kn 
iV 
2p7r 



pn 



{k + p)tt 



2{k+p)Tr, 
N ' 



(53) 



With these expressions, we are ready to compute the 
force correlation function in terms of dynamic variables 
Q^. We write 



F, 



k.fi 



N 



(54) 



P-2 



is the difference between bare 



where Aw^-^ 

dispersion relation and effective dispersion relation. The 
second term f^^ denotes the rest of the nonlinear force 
(we take only the quadratic terms in Q). Due to transla- 
tional and reflective symmetries, the correlation matrix 
formed by is diagonal without any approximation. 
The time-displaced correlation for the diagonal terms de- 
fines true force correlation. The nonlinear part of the 
contributions is 



E E 



ki-\-k2—k /c3+fc4 — 
YY YY^ 



cVm <^IIm (Qi (t)Qi (t)Qi* (^)Qi: (o)) > (ss) 
rT^..W«/3 E E 

c^k^. 4Il (qI ml: (^)Qi* (0)) • (56) 

In order to have a closed system of equations for the 
normalized correlation functions, we use the standard 
mean-field type approximation, (QQQQ) ~ {QQ){QQ)- 
Owing to the S correlation between different fc, the double 
summation can be reduced to a single one. We introduce 



<k{t) 



N 



kit)/p' 



277 A; 
'Na' 



(57) 
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FIG. 1: Damping functions of the slowest modes, 
part of rf [z] (solid line) and ri"[z] (dashed line) vs z, com- 
puted from molecular dynamics for parameters of set E 
Kr = 1, = 0.3, a = 2, T = 0.4, m = 1. The top fig- 
ure is for size N = 8 and bottom figure for N = 64. 



and similarly (i) associated with F^. In terms of z^a? (t), 



we obtain 



ki-\-k2—k 



ki+k2=k 



where 



K 



ki ,k2 



K 



ki.k-2 



2kBT 



knT 



27r(fci+fc2) .-._L - _L 
^kt^k2 

MY 
^ki,k2 



Na ki k2 



(58) 
(59) 

(60) 
(61) 



Equations (|^ 



and H59|) together with the relations 

among J^^f fc' ^ fe' ^'^'i -f^^^- G3' 63' ^'^^ (|S7jl . form 
a system of close equations, which can be solved in prin- 
ciple. However, because of the singular nature at z = wo 
in Eq. (|42|l . any approximation to Fjv will destroy a sub- 
tle cancellation of the singularity, rendering the problem 
impossible to solve. 

The damping function T'^lz] is the central function 
that a successful theory needs to be able to calculate. 
In Fig. ^ we present examples of such functions deter- 
mined from equilibrium molecular dynamics (MD) simu- 
lation in a microcanonical ensemble with periodic bound- 
ary condition. A more correct comparison of MD with 
mode-coupling theory should use an ensemble of ini- 
tial conditions distributed according to canonical weight. 
This may be unimportant when N is large. We com- 
pute F[z] from the ratio of two correlation functions. 



9Qq[^]/9Qp[^] = iz + T[z]. For small systems, there are a 
lot of peak structures associated with the low frequency 
modes; the feature appears washed out for large systems. 



V. SOLUTION OF MODE-COUPING THEORY 

A. Numerical solution at finite TV 

To make the problem tractable, we make a bold ap- 
proximation. Instead of Eq. H43|) . we take 



(62) 



This is equivalent to say Co ~ luq- This appears justi- 
fied for the longitudinal modes at sufficiently low tem- 
peratures but problematic for the transverse mode, as 
Lu is linear in wave number q but luq is quadratic in q. 
We can consider the limit of small q. In such limit, the 
difference between Tp and F is dropped. Lepri's treat- 
ment is similar to the above approximation. We also 
note that in the work of Scheipers and Schirmacher for 
damping of anharmonic crystals [s^ , effective cubic cou- 
pling is used instead of the 'bare' coupling, Ck,p,q- The 
standard operation of projecting the random force onto 
bilinear form QpQq to get the mode-coupling equations 
also agrees with the approximation, Eq. (|62|l . but with 
somewhat different "vertex functions" replacing Eq. H6(Jfl 
and 

In order to obtain a numerical solution, besides the 
model parameters (mass m, lattice spacing a, couplings 
Kr and K^), we also need the effective dispersion rela- 
tion. We used MD data for this purpose. It turns out 
that a two-parameter fit of the form 



2c I fcTT I 

ujk = — sm — 

a ' N ' 



- sm^- (63) 



characterizes the effective dispersion relation very well, 
where c is sound velocity at small k and Wmax is maximum 
frequency at fc = N/2. 

We used a fast Fourier transform to solve the equa- 
tions iteratively. Given an initial viz], correlation func- 
tion in frequency domain is calculated by Eq. In- 
verse transform gives us g{t), which is used to calculate 
the discrete sum in Eqs. (|58|l and H59|) to obtain i'{t). We 
transform back to frequency domain for the next itera- 
tion. 

Figure is our theoretical calculation with parame- 
ters exactly at set E. This figure should be compared 
with the MD data in Fig. ^ We observe that the pre- 
dicted results are too smooth without much structure. 
The mode-coupling theory predicts a damping that is a 
factor of two larger for the longitudinal modes and much 
too large for the transverse modes. By adjusting the 
temperature to a lower value (one half for N = 8 and 
a quarter for N = 64), we obtain curves which closely 
resemble the MD results. However, the mode-coupling 
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FIG. 2: Real part of Vl [z] (solid line) and Tjr [z] (dashed line) 
vs z from a full mode-coupling theory for parameters of set E 
(same as that in FigQ at A'' = 8 and 64. The input effective 
dispersion relation parameters, c" , c"'", ojmax, and Lu^a,^, are 
given in Table Q 
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FIG. 3: Real part of Tl [z] (soUd hne) and Fj^ [z] (dashed line) 
vs z from a full mode-coupling theory for parameters of set 
E, but at different temperature T. Other parameters are the 
same as in Fig. |21 



theory still overestimates the value of transverse damp- 
ing. It is clear that K-^ is considerably smaller in actually 
system, but our mode-coupling theory can not produce 
a small enough K-^. 



Set E bare 


TV = 8 


iV = 64 




2 


1.435 


1.341 







0.621 


0.669 


^max 


2 


1.495 


1.543 


^max 


1.095 


1.087 


1.182 


vo 


2 


0.7113 


0.6333 


Vl 





-0.1060 


-0.0221 


V2 


-0.6 


-0.1099 


-0.1242 


V3 





0.6362 


0.8965 


V4 





-0.0424 


-0.1346 



TABLE I: Parameters for mode-coupling equations. 

TV = 1024 

1.329 
0.674 
1.553 
1.194 
0.6337 
-0.0076 
-0.1293 
0.9312 
-0.1469 



B. An effective Hamiltonian approach 

The results presented in Fig. Ogive us some hints that 
the mode-coupling equations are essentially correct, but 
the parameters of the model need to be adjusted. In 
fact, simple perturbative expansion of the Hamiltonian 
can not correctly predict uj^. According to MD result, 
uj^ oc k, but the leading contribution from perturbation 



calculation should give the 'bare' frequency, tD^ = uj-^ oc 
P. 

The adjusting of the parameters can be made more 
rigorous (SQj with an introduction of an effective Hamil- 
tonian, 



fe,,.=|l,_L 



k+p+q=0 

- H" + H'. 



k+p+q=0 



(64) 



The form of the Hamiltonian is dictated by symmetry. 
Translational invariance requires that k + p + q — 
(mod TV). The system is symmetric under reflection 
about y-axis or x-axis. Thus the Hamiltonian should be 
invariant under the transformation '^Qk Q\ ~^ 

— qI^,. In addition, V^pq is symmetric under permutation 

~ (3) 

of p and g, and V^^^ is symmetric under permutation of 
all three indices. We also have V-k-p-q — —Vk^p^q and 
similarly for V^.^^. Although the original Hamiltonian 

does not have V"'^^ term, such a term can present. One 
of the reasons that such term can present is due to the 
non-analytic behavior of the potential (because of the 
absolute value | • |). 

We determine the parameters of the effective Hamil- 
tonian, (D)!, w^, V, and V^^\ by fitting them to the ob- 
served time- independent correlation functions from MD. 
In order to be able to carry out the calculation analyt- 
ically, we treat the interactions as perturbations, i.e., V 
and 

are assumed small. All the equilibrium aver- 
ages are approximated by the leading contribution from 
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the perturbation, 



0.20 



{l-f3H'))o 



(65) 



where (• • • )o is thermodynamical average with respect 
to the non-interacting harmonic oscillators (product of 
Gaussian integrals). Under the above approximation, we 
find 



HQ' 



M|2\ _ 



1 



and the interaction parameters 



kpq 



iQlQpQi) 



2P{\Ql\^){\Q^\^)(\Q^\ 



(3) 



{QlQlQl) 



kpq 



mQl?){\Ql?){\Q\ 



(66) 



(67) 



(68) 



The actual form of the u is fitted according to Eq. (|63|1 . 
For V and V'^'^^ wc fit into a functional form 



^kpq — 



ixyz 



{vo + wi + V2 {x^ + y^)) , 



~ (3) ^ ixyz 

kpq 



(u3+w4 (x^+y^ + z^)), 



(69) 
(70) 



where z = sin(fc7r/Af), x ~ sm[p'!i/N), and y — 
sin^qiT / N) . The values of vq to W4, together with c and 
<^max are listed in Tabled It is surprising that not only is 
the y*-^-* term present, but also is its magnitude compara- 
ble to V. The column marked 'bare' are the parameters 
corresponding to the 'bare' Hamiltonian, Eq. (|44|l . The 
bare parameters are reached only at very low tempera- 
tures, such as T = 0.002. Since we have factored out 
the leading size dependence, we expect the parameters 
weakly depending on size N . 

Before presenting our results with the effective param- 
eters, it is worth pointing out that the standard proce- 
dure of deriving mode-coupling equations gives identical 
result as our approximation, Eq. I|62|) . The standard pro- 
cedure is to project the random force, Rk — to^Qk + Qk, 
onto a bilinear form of the basic variable Qfc. Applying 
the general projection operator, Eq. (O, to our case, this 
bilinear projector is 



QiQj 



{\Q^C 



Y^N^QrQl, (71) 



where i or j is a pair of indices, e.g., j = (fc,/i), and j 
is (— fc,/i). In evaluating the four-point correlation hmc- 
tions, we have neglected the perturbation term. It turns 
out that the last term projects out only the fc = com- 
ponent. Thus, the specific values of Ni are not needed. 

After applying projector V2 onto i?^, we neglect the 
difference between normal dynamics e^^ and anomalous 
dynamics for R'^it), e*^ We have 



/3{{e'^V2Rm)('P2Kml- (72) 
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FIG. 4: Real part of Vl [z] (solid line) and F^ [z] (dashed line) 
vs z from a full mode-coupling theory for set E, using effective 
parameters given in Tabled 



The mode-coupling equations are then, 

p+q=k 

E K^A{t)9i{t), (74) 

p+q=k 



where 



K. 



p+qQ-pQ-q)\'^ 



pq 



2 m\'){mn 

J{Rp+qQ-pQ-q)\ 
{\QI\'){W) 

f3 \{Rl+qQ-pQ-q) 

2 (iQlipxigfp) 



2 

4 



Vp.-p-q.q 



18 

1 



^P^q 
p.(3) 

-p-q-.p-M 



,(75) 
, (76) 
.(77) 



UlpLOq 



Note that the vertex coupling (RkQpQq) = tjJ^iQkQpQq) 
is proportional to the three-point correlation indepen- 
dent of the specific form of random force. This version 
of mode-coupling equation need not have the limitation 
of small oscillations, and can be applied to high temper- 
atures as well. 

Figure 0] is a calculation of the real part of F^ [z] for 
k ^ 1. Comparing with Fig. |21 and |3| we see that using 
effective parameters brings into much better agreement 
with the MD data. In fact, comparing to Fig. ^ most 
of the features are reproduced, such as the locations and 
heights of the peaks. The most important improvement 
is the ratio of parallel to perpendicular damping. 

In Figure we compare the normalized dynamic cor- 
relation functions for size N = 256 and the slowest mode 
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FIG. 5: The normalized correlation functions (t) (upper 
part) and gi (t) (lower part) on a A'^ = 256 system for data set 
E. The solid curves are from full mode-coupling theory, while 
the dotted dashes are from equilibrium molecular dynamics. 



k = 1. The frequencies of oscillations are slightly dif- 
ferent in MD and mode-couphng calculation, thus there 
are phase shifts at long times. The longitudinal damping 
agrees with each other extremely well. A fit of the log- 
arithm of amplitudes (maxima and minima) versus time 
for the MD data gives the decay rate 7!' « 0.00064, « 

0.000023, while mode-coupling predicts 7^' ~ 0.00066 
and -f^ w 0.000078, respectively. Mode-coupling theory 
with effective parameters still overestimates the trans- 
verse damping by a factor of 2 to 3. 

In Fig. |H1 we show the damping of the mode in terms 
of the decay constant 7^, which is defined by the ex- 
ponential decay in g{t) « cos{ujkt)e~'^''* . The symbols 
are obtained from MD simulation, while the curves are 
from the full mode-coupling theory with effective param- 
eters. For MD data, 7fc is obtained by a least-square fit 
to the amplitudes. For mode-coupling theory, we use an 
approximation, 7^ = iffepfc]. The bar means we aver- 
age over a window centered around ujk with width about 
7fc. We observe that excellent agreement between MD 
data and mode-coupling theory is obtained for the paral- 
lel component. The mode-coupling theory overestimates 
the transverse mode damping by some constant factor. 
In any case, the slopes of the curves agree very well with 
the expected results (to be discussed in the next subsec- 
tion), 7[! cx fc^/^ and 7^ (x k^. 



C. Large A'^ limit 

To simplify the equations further, we consider the large 
size limit. Since the large size asymptotic behavior is the 
most interesting, such limiting results are justified. We 



FIG. 6: The decay rate 7^ for each mode k for system size 
= 1024. Circles (7JJ) and triangles (7*!") are from molecular 
dynamics by fitting the correlation function, while the contin- 
uous curves are from a full mode-coupling theory for set E, 
using effective parameters given in Table Q 



thus consider the limit of p and keep only leading 
contributions in lattice momentum p — 2'Kk/ [No). In 
this limit the three kernel functions K'^_^ become con- 
stants. In addition, we define a A:-independent v function 
by taking the limit p — > 0. The leading fc-dependence in 
Vk is k^ ■ When this term is factored out, we expect that 
i/fc becomes independent of k in the limit of large N and 
small p. With these simplification and approximation, 
and converting the discrete summation to integral, we 
have in the limit of large sizes: 



J{t) = 

where 
K 



2^ 
2ti 



■K / a 
7r/a 



(aHI 9^^(^)2 + ^(3)^(^)2)^(78) 

(79) 



— 7r/a 



dqg\{t)g^{t), 



2/3c-L4m3 



K^^^ = 0, 



(80) 



if we use the perturbation expansion Hamiltonian, or 



if" = 4 



fay _vl_ 

\2) /3cll2 



X(^) = 36 



7)2 
^3 



/3cll 



2„_L2 ' 



(81) 



if we use the effective Hamiltonian. We also linearize 
the dispersion relation so that lj^ ~ c^p, = c" or 
are effective sound velocities for the longitudinal and 
transverse modes. 
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We shall refer to Eqs. (|78|l and (|79|l as simple mode- 
coupling theory and the finite N version as the full mode- 
coupling theory. 



D. Asymptotic analytic solution 

Assuming that the simplified mode-coupling equations 
represent the essence of the physics regarding the damp- 
ing and time-dependent correlation functions, we now 
consider analytic solution to Eqs. H78I) and H79|) . First 
we notice some constraints on the functions. Because 
v'^it) is a real function, we must have ^[z]* = vl—z]. 
Since g^(0) = 1, we must have i'^(O) = K^^/a, for v = \\ 
or _L (assuming K^^'' — 0). This implies that the Fourier- 
Laplace transform must be integrable: 



1 f°° K 
z/(0) = - / i^[z]dz = — . 



(82) 



only look at the dominant term in time dependence of 
i'll(i). For 5^(i)^, since large g-mode decays rapidly, we 
ignore the term with amplitude proportional to g^, and 
also drop the oscillatory term and approximate, .g^(i)^ « 

ie~^9 After integrating over q (we also extend the limit 
as from — oo to oo), we obtain the leading t dependence, 
as 



II K» 1 

4 \ TTi^n i 



(86) 



If the oscillatory terms are kept, they only contribute an 
exponentially small term, e"*^^ ■* ; and the terms 
give a contribution that decays much faster, as t~^^^. 
The contribution from JsT^^^ term can be neglected, be- 
cause it decays as (due to Eq. H91|l 'l. But this term 
does provide a crossover from z^^/"^ for intermediate z 
to the asymptotic z^^/^ for very small z. The Fourier- 
Laplace transform gives 



Next, the large z behavior can be obtained by integrating 
by part few times, using the boundary conditions g{0) = 
1, g'{0) = 0, g"{0) = -0)2, and g{t ^ oo) = 0. For 
simplicity, we have omitted the lattice momentum index 
q and mode index || or _L, since it is true for any gf^. With 
these results, we find 




b{izy 



-1/2 



(87) 



We now get an expression for the longitudinal correla- 
tion function. Formally, it is given by the inverse trans- 
form 



it), 



^ J dq gl{t)gl{t)e-^^'dt 



K 



K 



. ... 3 !^'dq{i:^l+i^l) + 0{\),m 
laz ZmZ'^ ./-Tr/a ^ 



where for i^", g^ and g^ are all 5^, while for ^ it is (7II 
and g^ . 

We can establish that 



cx (iz)' 



-1/2 



V [z\ cx const. 



(84) 



.9ll(i) = 



1 

2^ 



-iz - Tq[z\ 



■^2 _ (,\\2q2 _ izTq[z\ 



dz. 



where Tqlz] 
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^z] = bqyVt^, b = \k\\I^. The 



dominant contribution is from z when the denominator 
is close to zero. The integral can be approximately es- 
timated by the residue theorem. By location the poles 



(89) 



- c"^g^ - ^izbq"- ^ 



we obtain the dispersion relation for g\(t\ In the limit 
of small g, we find 



in the limit of small . We proceed as follows: first we 
assume that v^\z\ = or v^(£) = 2vQb(i). We then 
derive an expression for v^ify and i'"'"(t), and show that 
indeed v^(t) approaches a (5-function in a proper limit. 



When 



is a constant, the inverse Fourier transform 



for the correlation function g\z\ can be performed exactly, 
to give the transverse correlation function as 




where 



kq^ ■, and uj = 



gj.2g2 



r^V4 



(85) 



c^q. 



We have linearized the dispersion relation, uj k, c q. We 



z w d^q + i-foq 



3/2 



where 70 = j^K^^ /yOhy^. Therefore, 



g^it) « e"'"'*'''*cos(cllqt). 



(90) 



(91) 



Similarly with the dispersion relation for the transverse 
mode, z w c^q + ^iv^q^, we can compute 



pOO 




cos{c" qt) coa{c qt) 



(92) 



We have neglected the 70 term. Since Ic" —c-^ \ is smaller 
than c" -\-c^ , we drop the second term. We symmetrically 
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extend the function. Note that i^'^{t) can be casted into 
the functional form 



2V7r ^ 
which has the property that 



(cll - c^f 



/OO 
f„{t)dt = l, for all <T>0. 
-OC 

Thus fcr{t) behaves as a ^-function as a — > oo: 



(93) 



(94) 



(95) 



This allows us to identify the constant = }^K^ l\c} — 
c^|. The condition for cr —> cx) is the same as Vq ^ Q 
or 0. Since K^jK^ = 2{c^lc}\ f (c.f. Eq.llHOl), 

small also implies small . 

Although the above derivation suggests that the result 
is valid only for the special limit. But in fact, this is 
also the asymptotic result in the limit z ^ 0. To show 
this, we note that the asymptotic behavior is picked up 
by the scaling lim^^o \^v\\z\ k, v[z\. If this is true, then 
y[z\ cx . This scaling limit coincides with the limit of 
0. 



E. Numerical solution of the simple theory 

For the simple mode-coupling theory that we have al- 
ready taken the large N limit, we can also apply the 
fast Fourier transform method. However, we find that a 
direct numerical integration in frequency space is much 
more accurate. For example, [z] is expressed as 



p+oo p+oo 

(p- + - 1. - 1.')) , (96) 



where P stands for Cauchy principal value and 



.9"'^('Z,^) = 



TT / a 



(97) 



(98) 



When the linear approximation is made to the dispersion 
relations, the integral over q can be performed analyti- 
cally. We only need to do a two-dimensional integral 
over w and w'. The principal value integral is taken care 
by locating exactly the singularity. Since the integration 
routine needs a smooth function v\z\, we fit the results 
by a Fade approximation (in variable iz or (iz)^/^). The 
procedure is iterated several times for convergence. This 
is programmed in Mathematica. It turns out essential to 



10' 
10' 



10 



10"' 



10" 



10 



10" 



10"' 




10' 



FIG. 7: Real part of ^/"[z] (top) and u^[z] (bottom) for 
data set E. The smooth curves are from simple mode-coupling 
theory, the (red) curves with spikes are computed from full 
theory at A'^ = 1024, while the straight, dash lines are analytic 
asymptotic results. 
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FIG. 8: Real part of iy"[z] (solid line) and u [z] (dashed line) 
for data set J (if" = 1.532, = 0.719, and A"*^' = 2.776). 
The insert shows the effective exponent ~d\nv /d\nz. 



have more than double precision accuracy in the integra- 
tion routine in order for the singular integrals properly 
converged. Some results of this calculation are already 
presented in ref. [3ll |. 

In Fig. we compare three levels of approximations of 
the simple mode-coupling theory, the full theory 



F^'l 



z]l{2^l{Na)y 

U _ 



computed on iV = 1024 by v^lz] 
and the asymptotic result of v'^[z\ — ^K-^il/\c'^ 
l/|cll-|-c"'"|), and Eq. H87|) . Noticeable differences are seen 
between full theory and simplified. This is partly due to 
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ing conditions 



FIG. 9: Real part (top) and imaginary part (bottom) of 
X^/'-^u{Xz; K, A^/^c, X-^^^a) vs z for if = c = a = 1 and A = 1 
(solid line), 1/8 (open circles), and 8 (crosses). 



the fact that we are comparing a finite N with an infinite 
N result. We also note that the asymptotic slope of —1/2 
for i^ll [z] is approached rather slowly. 

Figure |S| shows a stronger crossover effect from slope 
of 1/3 (corresponding to k (x N^^^, to be discussed later) 
to that of the true asymptotic law of 1/2. This set of 
curves corresponds to data set J. 



F. Scaling solution of the simple theory 

We now look at some of the scaling properties that are 
implied from Eqs. H78I) and (|79|l . First, we consider the 
simple case of if = if" = if-"-, K^^'^ = 0, and c = c" = 
c"*-. In this case the two equations degenerate into one 
equation given by Lepri 0. By measuring frequency in 
terms of c (i.e., consider variable z/c) we can scale away 
c, thus the following equation is an exact scaling 



v{z; K, c,a) — c i^{z/c; K/c^, 1, a), 



(99) 



where is a function of z with parameters if, c, and a 
that we have written out explicitly. This equation tells us 
that out of the three parameters that characterize the so- 
lution, we can pick two as independent. The "shapes" of 
the solutions are all the same for if/c^ = const. Without 
loss of generality, we can fix if = 1. 

Next, we consider a general scaling solution of the form 



i^(Az; if , A'^'^c, A'^°a) « A "2^(2; if, c, a). 



(100) 



Substituting this scaling form into Eq. (|95|l , by requiring 
that the result must be consistent, we obtain the foUow- 



A, 



1- A 



Aa + A„ = 



0, 
0, 
S, 
0, 



(101) 
(102) 
(103) 
(104) 



where Aq is scaling exponent associated with integration 
variable q, q ^ X^''q. A unique solution to the set of 
linear equations is obtained, S = 1/3, Ac = 1/3, Aa — 
—2/3, and A^ = 2/3. Since Eq. (|100|1 is (approximately) 
valid for any A, we can choose A = 1/z. This scaling 
solution implies that 

i^iz; if, c, a) « z~^/^iy{l; K, c/z^^^, z^l^a). (105) 

Power-law behavior z^^l'^ is obtained in the limit of small 
z, relatively large c, and small a, if i^(l, if, cxd, 0) is finite. 
The crossover to other behavior occurs at large z ^ c" 
and z ~ a""^/^. 

A simple dimension analysis also leads to the z^^l'^ 
factor. Let the dimension of length and time be L and T, 
respectively. Then the dimensions of relevant quantities 
are [z] = T^^, [c] = LT"!, [a] = L, {v\ = l?T'^ , and 
[if] = L^T"^. From the five quantities, we can construct 
three dimensionless variables 



Hi 



j-l/3^2/3' 



TI2 = — , Ha = (106) 
z a ac^ 



If there is any relation between the n^'s, it must be in 
the form Hi = /(n2,n3) (Buckingham Pi theorem 40J), 
or 



-1/3^2/3 J 



za if 



(107) 



where / is an arbitrary, dimensionless function. This 
result, of course, is consistent with Eq. H100|l . This also 
suggests that the scaling is not approximate, but an exact 
result. 

Fig. 1^1 is a test of the above scaling by numerical solu- 
tions. For small z, the power-law z~^/'^ is verified to high 
accuracy. Since this scaling is exact, the deviations are 
due purely to numerical errors in solving the equations. 

The case of the coupled longitudinal and transverse 
equations is somewhat difficult to analyze. We can re- 
quire a very general scaling of the form 

J (Az; X^kW if II , A'^'--^ if-L, A'^cii c" , A"^-^^ , A-^-a) 

« A"''" v{z; if II , K^,c} , c-L, a), (108) 

and similarly for the perpendicular component. If we 
require for consistency of scaling for both longitudinal 
and transverse components, we find that the only scaling 
solution is the symmetric solution, i.e., the scaling dis- 
cussed above with identical longitudinal and transverse 
scaling exponents. 
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If we abandon the exact scaling for the transverse com- 
ponent and look only at longitudinal component, we can 
require that 



Si_ = 2A„-1, 



6\\ = 1-A 
A,.± = 1-A„ 



(109) 
(110) 
(111) 



As before we can require that there is no scaling for 
the coupling constant i^" without loss of generality (i.e., 
A^ii — 0), the above equations imply a relation 



2Sn 



1. 



(112) 



In particular, if d± ~ 0, we must have 8\\ = 1/2 and 
A^± = 1/2. This is consistent with the fact that the 
scaling region of 5\\ = 1/2 is obtained for small . 

A complete and clear scaling picture is still lacking. 
From the numerical solution, we observe that besides well 
characterized scaling regimes, there are also plenty of in- 
termediate regions and crossovers. For very large and 
small K^, or small c?- and large K^, the behavior of the 
solutions are difficult to characterize. 



rent per particle: 

m], = - Ari((pj + Pi+i) • G(i)) 

-Ar,_i((p,+p,_i) -0(1-1)) 
+ Ar,_i(p, •H(i-2,i-l,i-l)) 
-HAri(p, •H(i-|-l,i + l,i)) +pi/ii, (115) 

where G(z) = \Kr{\lS,r^\ - a)nj, H(i, j, /c) = K^{n^ + 
rife cos 0j)/| Ar^l, is a unit vector in direction of 
Ari = r^+i — Yi. The total heat current in x direction, 
J — ^ j jx^i is the quantity appearing in the Green-Kubo 
formula. It is equal to the macroscopic heat current den- 
sity (energy per unit area per unit time) integrated over 
a volume. 

For theoretical analysis, we need to express J in terms 
of the dynamical variables Pk and Qk- A general ex- 
pression will be rather complicated. Again as in mode- 
coupling approach, we'll consider small oscillation expan- 
sions, and consider only the leading contributions. Ne- 
glecting the nonlinear contribution of ©(Qf), we obtain 



'dim 



(116) 



VI. GREEN-KUBO FORMULA 

In this section, we make the connection of the damping 
of modes with thermal conduction. The starting point is 
the Green-Kubo formula for transport coefficient: 



where is the bare dispersion relation given by 
Eqs. (|45|l and 1)46(1 . More general expression in 
a quantum-mechanical framework including the cubic 
terms is given in ref . [42j . 

An approximation for the correlation function of the 
current can be obtained, again using a dynamic mean- 
field approximation, as 



1 



{J{t)Jmdt. (113) (j(t)j(o)) =Y.\bt\'m{t)Q^_,m{Pk'it)PMO)) 



For the special case of heat conduction in one- 
dimensional chain, it is better to call k heat conductance 
rather than heat conductivity. In analogous to electric 
circuit, the heat conductance relates the temperature 
gradient to energy current (Fourier law). 



I=-K — 

dx 



(114) 



The quantity J{t) is related to energy current by J = 
laN. 

The central quantity that we compute in equilibrium 
and nonequilibrium heat conduction is the (total) heat 
current J. Since the heat current is a macroscopic con- 
cept determined by the conservation of (internal) energy, 
a microscopic version of it is not unique. The current 
expression must satisfy the energy continuity equation 
in the long- wave limit. We derive an expression start- 
ing from J = '^_^d{rihi)/dt, where the local energy per 
particle is hi = iif^ ((| Ar^^i | - a)^ + (|Ari| - a)^) + 
Kcj, cos +pf/{2m). By regrouping some of the terms 
using translational invariance, we arrive at the heat cur- 



k,^ 



+{Q'kit)PMo)r}. 



(117) 



The above expression can be further simplified by the 
approximation Pj^ = ~ lu'^Q'^. We find 



(J(OJ(0)) 



E 

k,fj, 



/3-fe" 



aiV(cll)^ 
27r/32 



g'kitf 

°°cos2(c:;k) e-^T^dp 



^i-i/(2-.„) 



(118) 



We have dropped the contribution from the perpendicu- 
lar component, because it decays much faster (due to an 
extra factor in the integrand) . 

We compare the mode-coupling result of Green-Kubo 
integrand with that of MD in Fig. ^| We have used 
Eq. (|117|l for the mode-coupling calculation. Note that 
the correlation functions {Qit)Q{0)*), (P(i)P(O)*), and 
(Q(t)P(O)*) are simply related in frequency domain 
through Eq. (|32|l . The agreement is reasonable with the 
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FIG. 10: The Green-Kubo integrand from MD (circles) and 
mode-coupling theory (curve) for set E with size A'^ = 1024. 
The straight has a slope —2/3. The insert shows the same 
data on a linear scale. 



(119) 



Since (Jy = 1/2, koc v.^- This result is in agreement with 
that of Deutsch and Narayan for a model (the random 
collision model) where transverse motion is taken into 
account only stochastically [iGj . 

In Fig. 111! the mode-coupling result of kat is compared 
with equilibrium molecular dynamics result with periodic 
boundary conditions. Excellent 1/2 power is observed for 
the mode-coupling result. We found it is rather difficult 
to get converged value from molecular dynamics. But 
the results are consistent with a \/N law. 

The relation of the mode-coupling theory with the re- 
sults of noncquilibrium situation of low and high tem- 
perature heat-baths at the ends is not clear cut. The 
standard assumption is that the correlation should be 
cut off by a time scale of order N due to the interaction 
with the heat baths. Finite size result is obtained by in- 
tegrating the power-law decay to a time of 0{N). This 
gives us the asymptotic behavior for the conductance at 
large N as 



100 



10 




10 



100 



1000 



(120) 



Since we find i5|| =1/2 for small , the thermal con- 
duction diverges as N^/'^. When K'^'^^ is sufficiently large, 
we observe N^^^. 



VII. NONEQUILIBRIUM MOLECULAR 
DYNAMICS STUDY 



A. MD simulation details 

We note that the interaction potential is not smooth 
at the point when two particles overlap. This can cause 
numerical instability. Thus, we replaced the original po- 
tential with a modified one, 



FIG. 11: The finite lattice conductance defined by Green- 
Kubo formula at parameter set E. Solid circles are from full 
mode-coupling theory; the triangles are from molecular dy- 
namics with periodic boundary condition; the straight lines 
have slope of 1/2. 



largest deviation about a factor of 2. The asymptotic 
behavior of t^^/-^ is consistent with both sets of results. 

The heat conductance on a finite lattice is sensitive 
to boundary conditions. Clearly, in the mode-coupling 
formulation, we have used periodic boundary condition. 
If we integrate over time from to cxd first on the second 
line of Eq. (|118|) . we obtain k cx Jdp/jp. On a finite 
lattice we have a lower momentum cut-off, 2tt / {No). The 
size dependence of the conductance on finite lattice is 



Ar = |rj+i - r.; 



Ar' = At 



Ar 



(121) 



which smooths out the discontinuous derivatives at Ar = 
with a small correction of order e^. In addition to 
replacing the spring potential by ^/^^(Ar' — of', we also 
replace the cosine term by 



COSffl,, = — ■ 



Ar,_i ■ Arj 
Ar^.iAr^ 



(122) 



so that the value is well-defined for all positions r^. 

In actual simulation, we have used very small e ^ 10"'^ 
to 10"'' so that its effect should be comparable to error 
caused by finite time-step h. We also used large e as a 
way of simulating a slightly different model to study the 
robustness of the results. 
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We solve the system of equations |4l| 



0.50 



dt 




if i < 

if i > N- N,„- 



(123) 



where is the total force acting on the i-th particle, 
and £,H obey the equation 



di 



L.H 



1 



dt 



kBTL,HN„ 



y El 



(124) 



where and are the temperatures of the two heat 
baths at the ends. The summation is over the particles 
belonging to the heat bath. We have used 4 particles for 
the Nose-Hoover heat baths, with extra first two and last 
two particles at fixed positions. The coupling parameter 
is taking to be 1. For the central part, we used a 
second-order symplectic algorithm (or equivalently, the 
velocity Verlet algorithm), while for the heat-bath, we 
used a simple difference scheme accurate to second order 
in h for ^. 

Although the total energy (Hamiltonian) is no longer a 
conserved quantity when the heat-bath is introduced, the 
above equation still has a conserved quantity of similar 
character: 



H{p,q) 



E 

•=L,H 



^xdt . (125) 



This quantity can be used to monitor the stability of the 
algorithm. 

Since we run for very long time steps of 10^ to lO^*', 
it is essential that the algorithm is stable over extended 
period of time. Even with a symplectic algorithm, stabil- 
ity is not guaranteed by merely taking small h (^ 10^"*). 
Thus, we adjusted the conserved quantity, Eq. (|125|l . to 
its starting value after certain number of steps. 
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FIG. 12: Temperature of i-th particle, Ti vs scaled position, 
i/N , for data set E with N = QA (plus), 256 (dash line), and 
1024 (solid line). 
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B. Heat conductance results 

In Fig. 1121 we observe that good temperature profiles 
are established. Due to relatively large temperature dif- 
ference between low and high temperature heat baths, 
and perhaps also due to the nature of heat baths, the 
profile is not linear. However, the scaling with N is ap- 
proximately obeyed. 

In addition to the thermal conductance results re- 
ported in ref. ,31,], Fig.ll3lgives additional data for variety 
of parameters. One of the aim of these additional runs is 
to check the robustness of the 1/3 power law for thermal 
conductance. It is found in ref ^21] for set E exponent 
a — 0.334 ± 0.003. This is an excellent confirmation of 
the mode-coupling theory. We note that at parameter 
set K, the 1/3 law is not destroyed by introducing large 
e. Thus we believe that the 1/3 law is not specific to the 
potential used in ref. [IJ. In fact, the Fermi-Pasta-Ulam 



FIG. 13: jN vs iV for different parameters {K4,, Tl,Th, e). set 
C: {K4,,TL,TH,e) = (0.1,0.2,0.4,0), the chain is compressed 
to have an average distance between particle do ~ 1.5; set G: 
(10, 0.2, 0.4, 0); set I: (0.1, 0.3, 0.5, 0.2); set K: (0.5, 1.2, 2, 
0.4); set L: (25, 1, 1.5, 0.2). All of the sets have Kr = 1, mass 
m — 1, and lattice constant a = 2. The number indicates the 
slope of least-squares fit. 



model with transverse motion was studied [43 and result 
is also consistent with the 1/3 power law. 

At low temperatures for set I and C, we appear to ob- 
serve exponent of 0.4. What is particularly interesting 
for set C is that the chain is compressed to an aver- 
age distance of 1.5 from the equilibrium distance of 2. 
The behavior is not changed much comparing to uncom- 
pressed chain. Set G and L represent very large angular 
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coupling K(f,. We note that K^j, oo corresponds to a 
situation that equihbrium can not be estabhshed. As 
becomes larger, the equilibration becomes increasingly 
difficult. In ref. '3l] we reported a logarithmic behavior 
for Kef, = 1 (set B). This logarithmic behavior does not 
maintain when K(j, is increased further. The exponent a 
for set L and G is closed to 1/2. This may be related the 
results of disordered harmonic chain. 



VIII. CONCLUSION 

Lepri et al first introduced mode-coupling theory into 
the problem of heat conduction in low dimensional sys- 
tems to interpret the MD results qualitatively. We have 
developed the theory further by introducing a 'full the- 
ory' and effective parameters for theory. The full theory 
is solved numerically. At this level of approximation, we 
find that mode-coupling theory gives excellent prediction 
on a quantitative level. The longitudinal damping (de- 
cay rate) agrees with MD data with few per cents of 
deviation. The mode-coupling theory somewhat overes- 
timated the transverse damping 7^ by a factor of two to 
three. The full theory assumes that the damping func- 
tion (memory kernel) Tk [z] is an explicit function of two 
variables, k and z. In a simple mode-coupling theory, we 



assume Tk[z] « {2Trk/{Na))'^h'[z], valid for small k/N. 
The simple theory is computationally efficient, and con- 
tains essential asymptotic features, i.e., j^, oc k^^'^ and 
7^ cx /c^. Through detailed comparison between MD and 
mode-coupling theory on memory kernel, decay of the 
modes, the time-dependent correlation functions, Green- 
Kubo integrand, and finally the finite-size conductance, 
we have strong evidence that the mode-coupling theory 
is correct, and it captures the essential features of the 
system. 

We have used parameter set E as the main data set 
for comparison. The reason for choosing this particular 
set of parameters is that it gives the cleanest power-law 
behavior. Other parameters will be either qualitatively 
similar, or crossovers will be observed. In fact mode- 
coupling theory naturally predicts crossover, due to the 
presence of K^^^ term. When this term is dominant, 
or the transverse effect can be neglected, we should see 
behavior of k oc 7V2/5. Thus we feel that the controversy 
of the result of Lepri et al m and that of Narayan and 
Ramaswamy (k cx N^^^) |15| can be reconciled within 
the current mode-coupling theory. 
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